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ABSTRACT 

This  report  describes  the  method  of  integration  that  is  used  in  a  Lincoln  Laboratory 
computer  program  (the  Planetary  Ephemeris  Program)  to  determine  as  functions  of  time 
the  position  and  velocity  of  the  Moon  and  the  partial  derivatives  of  these  quantities 
v/ith  respect  to  initial  conditions.  The  method  consists  of  numerically  integrating  the 
differential  equations  for  the  differences  between  the  positions,  velocities  and  partial 
derivatives  in  the  true  lunar  orbit  and  in  Brown's  mean  lunar  orbit. 
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GENERATION  OF  THE  LUNAR  EPHEMERIS 
ON  AN  ELECTRONIC  COMPUTER 


I.  INTRODUCTION 

A  computer  program,  called  the  Planetary  Ephemeris  Program  (PEP),  is  being  written  at 
Lincoln  Laboratory.  The  purpose  of  the  program  is  to  improve  planetary  and  lunaj-  ephemerides 
by  using  the  results  of  radar  and  optical  observations.  In  this  report  we  describe  the  method  of 
integration  that  is  used  in  PEP  to  determine  as  functions  of  time  the  position  and  velocity  of  the 
Moon  and  the  partial  derivatives  of  these  quantities  with  respect  to  initial  conditions. 

In  Ref.  1  we  described  Encke's  method  of  integration  used  in  PEP  for  the  planets  and  the 
Earth-Moon  barycenter;  we  also  presented  Encke's  method  of  integration  for  the  Moon.  How¬ 
ever,  this  method  is  not  satisfactory  for  the  Moon,  since  the  motion  of  the  Moon  deviates  greatly 

from  elliptic  motion.  In  fact,  the  expression  given  in  Brown's  lunar  theory  for  the  mean  lunar 
2 

orbit  implies  that  the  Moon  approximately  follows  an  ellipse  whose  ascending  node  moves  back¬ 
ward  along  the  ecliptic  one  revolution  in  18.6  years  and  whose  perigee  advances  one  revolution 
in  6  years. 

Our  method  of  integration  for  the  Moon  utilizes  Brown's  mean  lunar  orbit  rather  than  the 

1  A 

initial  osculating  elliptic  orbit  of  Encke's  method  of  integration.  Namely,  let  (x  ,  . .  . ,  x  )  denote 
the  position  and  velocity  in  the  true  orbit  of  the  Moon  and  let  (y1,  . .  .  ,y^)  denote  the  position  and 

velocity  in  Brown's  mean  lunar  orbit.  We  then  numerically  integrate  the  differential  equations 

1c  k  k  k  * 

for  the  |  =  x  —  y  ,  determining  the  x  from  the  results  of  the  integration  by  the  fact  that  the 

yk  are  known  as  functions  of  time.  The  partial  derivatives  dx^/dp^  of  position  and  velocity  with 

respect  to  initial  osculating  elliptic  orbital  elements  (/?*, . . . ,  p^)  needed  in  fitting  to  observations 

are  determined  either  (1)  by  assuming  that  these  quantities  are  equal  to  the  partial  derivatives 

3y  /dp  of  position  and  velocity  in  Brown's  mean  lunar  orbit  with  respect  to  initial  mean  orbital 
_ 1 

elements  ( p  ,...,/?  )  or  [ '.)  by  numerically  integrating  the  differential  equations  satisfied  by  the 

quantities  tj  ^  =  3 x^/pp^  —  dy^/djfi.  Opti  a  (1)  might  gh  e  results  of  the  accuracy  required  by  the 

least-squares  process  of  fitting  to  observations.  Not  having  to  follow  the  exact  procedure  of 

k  * 

(2)  would  save  a  great  deal  of  computer  time.  The  partial  derivatives  ox  /3«  of  position  and 
velocity  with  respect  to  parameters  a  that  are  not  initial  conditions  (such  as  the  mass  of  the 
Moon  or  the  combined  mass  of  the  Earth  and  Moon)  are  determined  by  numerically  integrating 
the  differential  equations  for  these  quantities;  since  these  equations  are  given  in  Ref.  3,  we  do 
not  derive  them  in  this  report. 

We  intend  to  integrate  the  equations  in  _  coordinate  system  referred  to  the  mean  equinox 
and  equator  of  1950.0.  The  use  of  this  coord  ..ate  system  dictates  some  of  the  manipulations 
we  perform  in  Secs.  Ill  and  IV  on  Brown's  mean  lunar  orbit,  since  the  reference  angles  for  this 
orbit  are  given  relative  to  the  mean  equinox  and  ecliptic  of  date. 
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The  differential  equations  of  motion  of  a  given  body  about  a  central  bod}  can  be  numerically 
integrated  with  arbitrarily  given  accuracy  over  an  arbitrarily  given  period  of  time  if  a  small 
enough  interval  of  integration  is  used  and  if  enough  figures  can  be  handled  in  the  computations 
to  prevent  significant  accumulation  of  round-off  errors.  By  subtracting  a  mean  orbit  from  the 
true  orbit  as  we  do  in  the  case  of  the  Moon,  the  requirements  on  the  size  of  the  interval  of  in¬ 
tegration  and  on  the  number  of  figures  needed  are  less  stringent  than  if  we  just  worked  with  the 
original  equations  of  motion. 

The  electronic  computer  (an  IBM  360  Model  67)  we  intend  to  use  for  the  numerical  integra¬ 
tion  has  hardware  floating  point  arithmetic  operations  with  execution  times  of  a  few  microseconds 
which  handle  7,  16  or  32  decimal  places.  As  presently  envisioned,  PEP  will  make  certain  crucial 
computations  (such  as  taking  a  numerical  integration  step  or  calculating  mean  lunar  orbit  quan¬ 
tities)  with  32  decimal  place  accuracy.  Other  computations  (such  as  the  determination  of  the 
positions  of  perturbing  planets)  will  be  done  with  16  decimal  place  accuracy.  If  necessary,  the 
design  of  PEP  could  be  altered  so  that  crucial  computations  could  be  carried  out  to  64  (or  even 
more)  decimal  place  accuracy.  Such  higher  precision  floating  point  operations  would  have  to  be 
programmed  operations  rather  than  machine  operations  and  would  thus  require  a  great  deal  more 
computer  time. 

In  any  event,  PEP  can  be  designed  to  handle  enough  figures  to  prevent  the  significant  accumu¬ 
lation  of  round-off  errors.  The  only  question  is  whether  so  small  an  interval  of  integration  is 
needed  that  excessive  computer  time  is  used  in  integrating  the  motion  of  the  Moon  for  centuries 
with  the  accuracy  required  by  observations.  This  point  can  only  be  determined  by  computer 
experimentation.  Because  of  the  rapidity  with  which  improvements  are  being  made  in  electronic 
computers,  such  a  calculation  could  very  well  be  handled  within  a  decade  if  not  at  the  present 
time. 

Since  PEP  is  written  in  the  Fortran  IV  language,  only  slight  modification  of  the  program  will 
be  necessary  for  its  use  in  future  computers.  It  is  quite  easy  to  insert  in  PEP  the  effect  of  ad¬ 
ditional  forces,  since  it  is  the  logic  concerned  with  making  a  numerical  integration  step  and  ma¬ 
nipulating  input  and  output  which  makes  the  program  intricate,  not  the  specific  terms  on  the  right- 
hand  sides  of  the  equations. 

II.  EQUATIONS  OF  MOTION  AND  EQUATIONS  FOR  PARTIAL  DERIVATIVES 
WITH  RESPECT  TO  INITIAL  CONDITIONS 

We  make  the  following  definitions  concerning  subscripts: 
s  =  Sun 
e  =  Earth 
m  =  Moon 

c  =  Earth-Moon  barycenter  (center  of  mass  of  Earth-Moon  system) 
j  =  jth  planet  (j  =  1,  2,  4,  ....  9) 

12  3 

Let  y  denote  the  gravitational  constant  and  suppose  that  (x  ,  x  ,  x  )  is  an  inertial  coordinate 
system.  We  make  the  following  notational  conventions: 

k  th 

xg  =  k  coordinate  of  s,  etc. 
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k  k  k  th 

x.  =  x.  —  x„  =  k  coordinate  of  j  relative  to  s  so  that 
js  3  s  J 

k  k  , 

-  -xsj-  etc- 

r  .  =  r-  =  distance  between  s  and  j,  etc. 
sj  3s 

Mg  =  mass  of  s,  etc. 

M  =  M  +  M  =  mass  of  Earth-Moon  barycenter 
c  e  m 

k  k 

Let  Fg  and  Fm  denote  the  components  of  force  on  the  Earth  and  Moon,  respectively,  in  addition 
to  those  due  to  the  planetary  and  solar  attractions.  Then  by  Newton's  laws  of  motion  and  gravity 
we  have 


— f  =  tM 
dt 


:k  xk  xk 

+  yM  -^+yVM.Je+l_Fk 
m3  1  s  3  1  Li  i  3^  M  e 

r  r  J  r  g 

me  es  j  je 


£4 

dt2 


„„  xem  ,  „„  sm  ,  „„  m  ,  1  „k 

=  5 —  +  yM  — j —  +  y  / ,  M.  — L—  +  ^ —  F„ 

'  e  3  '  s  3  u  y  3  M  m 

r _  r _  .  J  r.  m 


k 


rjm 


k  =  1,  2,  3 


(1) 


me  ms  3 

Subtracting  the  first  equation  of  (1)  from  the  second  equation  of  (1),  we  obtain 
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d2xk 


(sr) 


s'  r 


me  ,  _k  .  ,,,k  .  (  \  „k  1  „k 

__  +  B  +y  +  Fm__-F 

me  m  e 


)  ,  k=l,2, 


(2) 


where 


Bk  =  yM 


es 


es 


k  , 

ms 

3 

msy 


.  _  IVI .  / 

j  sV 


M.  /x.1'  x.  ' 
jm  _  je 

T"  T 

vrjm  V 


k  =  1,  2,  3 


(3) 


k  k  /  k  / 

Let  H  denote  that  part  of  the  (F  /M  -  F  /M  )  term  due  to  the  higher  harmonics  in  the  grav- 
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itational  potentials  of  the  Earth  and  Moon,  and  hereafter  let  F  *  and  Fm  denote  the  components 
of  force  due  to  effects  other  than  these  higher  harmonics  and  the  planetary  and  solar  attractions. 
Then  (2)  can  be  put  in  the  form 


dxme  _  k+3 

dt  -  me 


dx 


k+3 

me 


=  ->*1s (sr)  ^T-  *  B“ +  *k  +  »k  +  (if-  Fm-ifFe) 

'  s  r _  '  m  e  ’ 


dt 


xk  =xk 
me  ome 


me 


xk+3  =  xk+3  when  t  =  t 
me  ome  o 


k  =  1,2,3 


(4) 
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i 


which  is  exactly  as  written  in  Ref.  4.  The  expression  given  for  H  in  Ref.  5  includes  the  effect 

of  the  second  and  third  harmonics  of  the  Earth  and  of  the  second  harmonic  of  the  Moon. 

Let  (/3 1 , . . . ,  /3  ^ )  denote  the  osculating  elliptic  orbital  elements  at  time  t  of  the  orbit  of  the 
m 

Moon  about  the  Earth.  Differentiating  (4),  we  obtain0 


dOxL/^i)  «*£!? 
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L  (JL 

3  K 


F  k  - 
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dp 
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when  t  =  t 


3/3 


m 


3/3 


m 


k  =  i,  2,  3 
j  =  1, ....  6 


Jr  Jr  Jr 

where  the  partial  derivatives  of  B  ,  'k  ,  H  are  given  in  Ref.  7. 
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HI.  MOTION  IN  BROWN'S  MEAN  LUNAR  ORBIT 

In  the  following  the  unit  of  time  t  is  measured  in  ephemeris  days.  Let 

i  =  inclination  of  Brown's  mean  lunar  orbital  plane  on  the  mean  ecliptic 
of  date 

ft  =  ascending  node  of  Brown's  mean  lunar  orbital  plane  on  the  mean 
ecliptic  of  date  measured  from  the  mean  equinox  of  date  along  the 
ecliptic 

«  =  argument  of  perigee  of  Brown's  mean  lunar  orbit  measured  along 

the  orbital  plane  from  the  ascending  node  on  the  mean  ecliptic  of  date. 

2 

We  then  have 


sin -£=  0.044886967  (i*5!145) 

ft  =  259"?183275  -  0!0529539222  (t  -  t*) 

+  1!557  X  10"12  (t-t*)2  +  5?0  X  !0'20  (t  -  t*)3  ' 
w  =  T1  -  ft  =  75T14628 1  +  0°1643500025  (t  -  t*) 

-  9!296  X  10"12  (t  -  t*)2  -  3!1  X  10‘19  (t  -  t*)3 

where  t#  is  the  time  at  the  epoch  1900  January  0.3  =  J.  E.  D.  2415020.0. 
in  radians,  then 


(6) 


If  we  measure  angles 


4 


L 


51  =  4.52360151485  -  9.24220294234919  X  10-4  (t  -  t*) 

+  2.717477645355  X  10~U  (t  -  t*)2  +  8.72664625997  X  lO-22  (t  -  t*)3 
w  =  1.31155002408  +  2.868508295626071  X  10"3  (t  - 1*) 

-  1.6224580726539  X  10-13  (t  -  tj2 

-  5.4105206811824  X  10-21  (t  -  t*)3 
From  these  equations  it  easily  follows  that 


(7) 


^  =  -9.24220294234919  X  10-4  +  5.434955290710  X  10-14  (t  -  t*) 
+  2.617993877991  X  10-21  (t  -  t*)2 
~  =  2.868588295626071  X  10-3  -  3.2449161453173  X  10-13  (t  -  t,) 
-  1.62315620435472  X  10-20  (t  -  t*)2 


(8) 


d251 

dt2 


d2co 


dt 


=  5.434955290710  X  10-14  +  5.235987755982  X  10~21  (t  - t*) 


-3.2449161453178  X  10 


-13 


3.24631240870944  X  10 


-20 


(t 


(9) 
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Let  (v  ,  v  ,  v  )  be  a  coordinate  system  such  that  the  v  axis  points  toward  the  perigee  of  the 
2 

mean  lunar  orbit,  the  v  axis  lies  in  the  orbital  plane  and  points  in  the  direction  of  motion  at 

3 

perigee,  and  the  v  axis  is  perpendicular  to  the  orbital  plane  and  completes  the  right-hand 
12  3 

system.  Let  (w  ,  w  ,  w  )  be  a  coordinate  system  referred  to  the  mean  equinox  and  ecliptic  of 

1 

date,  that  is,  a  coordinate  system  such  that  the  w  axis  points  toward  the  mean  equinox  of  date, 

3  2 

the  w  axis  is  perpendicular  to  the  mean  ecliptic  of  date  and  points  to  the  north,  and  the  w  axis 

completes  the  right-hand  system.  The  relation  between  the  (v4,  v2,  v3)  and  (w4,  w2,  w3)  coor¬ 
dinate  systems,  assuming  that  they  have  the  same  origin,  is 


3  3 

wj=  £  B^vk  ,  vj=  l  Bkwk  ,  j  =  1,2,3  (10) 

k=l  k=l 


1  C 

where  the  orthogonal  matrix  B  =  (B^)  is  given  by 


B .  =  cos  51  cos  w  —  sin  51  sin  w  cos  i 
1 

1 

B£  =  —  cos  51  sin  w  —  sin  51  cos  w  cos  l 
1 

B^  =  sin  51  sin  i 
2 

B,  -  sin  51  cos  w  +  cos  51  s;nw  cos  i 
1 

2 

B^  =  —sin  51  sin  w  +  cos  51  cos  w  cos  i 
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Bg  =  —  cos  12  sini 


=  sinw  sini 


B 2  =  cos <v  sini 


B^  =  cosi 


(11) 


We  have 


dB 


dt 
d2B,:j 
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dt 
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912 


dev 
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(12) 
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(16) 

9  w 
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j,  k  =  1,  2,  3 


(17) 
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Let  (y  ,  y  ,y  )  be  a  coordinate  system  referred  to  the  mean  equinox  and  equator  of  1950.0, 

that  is,  a  coordinate  system  such  that  the  y1  axis  points  toward  the  mean  equinox  of  1950.0,  the 

y3  axis  is  perpendicular  to  the  mean  equator  (of  the  Earth)  of  1950.0,  and  the  y2  axis  completes 

the  right-hand  system.  Thq  notation  1950.0  denotes  the  instant  near  the  beginning  of  the  calendar 

year  1950  when  the  longitude  of  the  mean  Sun  was  18*1  40m  so  that  1950.0  is  J.  E.  D.  2433282.423  ^ 

12  3  12  3 

The  relation  between  the  (y  .y  ,y  )  and  (w  ,w  ,  w  )  coordinate  systems  is  given  by 


yJ  =  £  A* 
k=l 


w 


w 


■  2 

k=l 


j  =  1,2,3 


(18) 
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wnere  expressions  for  the  orthogonal  matrix  A  -  (A^)  and  its  derivatives  dA/dt  and  d2A/dt2  are 
given  as  functions  of  time  in  Appendix  B. 
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Combining  (10)  and  (18)  we  see  that  the  relation  between  the  (y  ,  y  ,y  )  coordinate  system 

12  3  1 

referred  to  the  mean  equinox  and  equator  of  1950.0  and  the  (v  ,  v  ,  v  )  coordinate  system  with  v 

3 

axis  pointed  toward  the  perigee  of  Brown's  mean  lunar  orbit  and  with  v  axis  perpendicular  to 
Brown's  mean  lunar  orbital  plane  pointed  toward  the  north  is 


3  3 

yJ  =  E  Cjjvk  ,  v^  =  E  Ckyk  .  j  =  l,  2, 3 

k=l  k=l 

where  the  orthogonal  matrix  C  =  AB  is  given  by 

3 

d-  E  A/Bk  >  j,  k  =  1.  2,  3  . 

£=1 


09) 


(20) 


The  derivatives  of  the  matrix  C  are 


dC  _  dA  R  A  d]3 
dt  "  dt  a  A  dt 


d2C  _  d2A  R  A  dA  dB  ,  .  d2B 
dt2  "  dt2  ‘  dt  dt  A  dt2 


(21) 


We  are  interested  in  a  body  moving  in  Brown's  mean  lunar  orbit  with  position  coordinates 
1  2 

(v  ,  v  ).  The  position,  velocity  and  acceleration  of  this  body  in  the  coordinate  system  referred 
to  the  mean  equinox  and  equator  of  1950.0  is  then 


y>-  l  cfc* 

k=l 


^  -  z  ci  *  z 


Tit-  '“'k  dt 

k=l 


2  dd 

k  vk 


dt 


k=l 


d  C 


,2  ]  1  .  ,2  k  2  dd  .  k 

d  r  _  y  rj  d  v  2  y  _ k  dv_  ,  _ 

"  U  ^k  ,.2  +  L  dt  dt  +  Li 

dt  k=l  dt  k=l  k=l  dt 


*d 


k  vk 


j  =  1.2.3 


(22) 


2 

The  mean  anomaly  L  at  time  t  in  Brown's  mean  lunar  orbit  is 
L  =  (J  -  T'  =  — 63!895392  +  13^0649924465  (t  -  U) 

r  6?889  X  10"12  (t  -  t*)2  +  2?99  X  10-19  (t  -  t*)3 


(23) 


If  we  measure  angles  in  radians,  then 

L  =  -1.1151849673  +  0.22802713493961401  (t  -  t*) 

1  1.202357321699  X  10-13  (t  -  t*)2  +  5.218534463463 


X  10"21  (t  -  t;.()3 


(24) 
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From  this  it  follows  that 


~  =  0.22802713493961401  +  2.404714643398  X  10“13  (t  -  tj 
at  * 

+  1.5655603390389  X  10"20  (t-t*)2  .  .  (25) 

2 

=  2.404714643398  X  10-13  +  3.1311206780778  X  10'20  (t  -  tj 
dt 

2 

The  eccentricity  e  of  Brown's  mean  lunar  orbit  is  constant: 

e  =  0.054900489  .  (26) 

The  semimajor  axis  a  of  Brown's  mean  lunar  orbit  is  also  constant.  By  inference  from  Ref.  2, 
we  can  assume  that 

a  =  60.2665  equatorial  earth  radii  .  (27) 

Let  u  be  the  eccentric  anomaly  at  time  t  in  Brown's  mean  lunar  orbit.  It  is  determined  by 
solving  Kepler's  equation 

L  =  u  —  e  sin u  (28) 

1  2 

by  iteration.  Then  the  radius  distance  p  and  the  components  of  position  (v  ,  v  )  in  Brown's 

10 

mean  lunar  orbit  at  time  t  are  given  by 
p  =  a(l  —  e  cos  u) 

v1  =  a(cos  u  -  e)  •  •  (29) 

2  /  2 
v  =  a  v  1  —  e  sin  u  . 

From  (28)  it  follows  that 

du  _  _ 1 _  dL 

dt  ~  (1  -  e  cos  u)  dt  ' 

so  that  by  (29) 

dv^  _  _  a  sin  u  dL 
dt  (1  —  e  cos  u)  dt 

0  , - -  (3D 

2  /  2 

dv  _  avl-e  cos  u  dL 

dt  "  (1  -  e  cos  u)  dt  . 


(1  -  e  cos  u)  | 

-sinu  — y- 
i  dt 

a*il  -  e2 

d2L 

(1  —  e  cos  u) 

COS  U  p 

dt 

d  L  (cos  u  -  e)  ,  dL  ,2 

2  .2  dt 

dt  (1  -  e  cos  u) 


dv  avl-e  I  dL  sin  u  ,  dL  ,2  | 

—T  =  (1  -e"cos  u)  C0S  U  TiT  “  - - : - 71  (  dFM 

dt  l  dt  (1  -  e  cos  u)  J  J 

The  position,  velocity  and  acceleration  in  the  coordinate  system  referred  to  the  mean  equinox 
and  equator  of  1950.0  is  then  given  by  (22). 


IV.  PARTIAL  DERIVATIVES  IN  BROWN'S  MEAN  LUNAR  ORBIT 


Let 

a  =  sernimajor  axis  at  time  t  of  Brown's  mean  lunar  orbit 

e  =  eccentricity  at  time  t  of  Brown's  mean  lunar  orbit 

i  =  inclination  at  time  tQ  of  Brown's  mean  lunar  orbit  to  the  mean 
equator  of  1950.0 

£2  =  ascending  node  at  time  tQ  of  Brown's  mean  lunar  orbit  on  the  mean 
equator  of  1950.0  measured  from  the  mean  equinox  of  1950.0 

w  =  argument  of  perigee  at  time  tQ  of  Brown's  mean  lunar  orbit 

measured  from  the  ascending  node  on  the  mean  equator  of  1950.0 

l  =  mean  anomaly  at  the  initial  time  t  in  Brown's  mean  lunar  orbit. 

In  this  section  we  shall  derive  the  expressions  for  the  partial  derivatives  of  the  position,  velocity 
and  acceleration  in  Brown's  mean  lunar  orbit  with  respect  to  these  quantities. 

We  of  course  have  a  =  a  and  e  =  e.  Let 


M 

“  =  <YMs)  jjf-  . 

s 


(33) 


Then  the  mean  motion  in  Brown's  mean  lunar  orbit  is 

1/2  -3/2 
n  =  p  '  a  ' 


(34) 


If  the  orbital  elements  of  Brown's  mean  lunar  orbit  were  not  functions  of  time,  we  would  have 
the  mean  anomaly  L  at  time  t  given  by 

L  =  I  +  n(t  -  tQ)  . 

By  this  equation  and  expression  (24),  the  mean  anomaly  L  at  time  t  in  Brown's  mean  lunar 
orbit  can  be  written  in  the  form 

L  =  7  +  n(t  —  t  )  +  L,(t  -  l  )  +  L,(t  -  t  J2  +  L,(t  -  t  )3 

O  *  O  Ci  O  J  u 


so  chat 


3L 

3a 

3L 

3e 


3L  , 

w  =  1 


3  n(t  "  V 


(35) 


Then  from  (28)  it  follows  that 


3u 

3  - 

3a  “ 

-  2  a(l  —  e  c 

3u 

sin  u 

3e 

(1  —  e  cos  u) 

3u 

1 

3£ 

(1  —  e  cos  u) 

(36) 
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"g=r. 


Differentiation  of  (29)  with  respect  to  a,  e  and  t  gives 

av1  .  3  nsinuU-g  ) 

W  =(cosu"e)+  2  “(l-e  cosu)  S 


I - 2  3  n  V 1  —  e  cos  u  (t  —  t  ) 

il  —e  sinu  — 4  - r: - ; - 

2  (1  —  e  cos  u) 


.  2 

a  sin  u 
(1  -  e  cos  u) 


ae  sin  u  a  v  1  —  e  sin  u  cos  u 

j - ^  ( i  —  e  cos  u) 

v  1  —  e4 


a  sinu 
(i  —  e  ccs  u) 


9v_  _  a  v  1  —  e  cos  u 
9f  ( i  -  e  cos  u) 


For  any  parameter  o’  that  is  independent  of  time  (such  as  the  initial  orbital  elements),  we  have 
|rr  =  -|=  Thus,  either  differentiating  Eqs.  (31)  and  (32)  with  respect  to  a,  e,  i  or  differ¬ 
entiating  Eqs.  (37),  (38)  and  (39)  with  respect  to  time,  we  obtain 


d_  (dv  \  _  3  n  sin  u _ 1 

dt  \  9a  /  ”2(1— e  cos  u)  (1  —  e  i 

d  /9v2\  _3  n  n/i  —  e2  cos  u  _ 

dt  \  9a  /  2  (1  —  e  cos  u)  (1 

d  / 9v*\  a  sin  u  dL  [ 

dt  [  9e  /  "  ,,  ,2  dt  COS' 

'  (1  —  e  cos  u)  «■ 

d_  / 9v2\  _  a  dL  _  e  co 

dt  \  9e  /  “  ( 1  —  e  cos  u)  dt  / — 


l _  dL 

cos  u)  dt 


3  n(t  —  t  )  (cos  u  —  e) 
sinu  — 2  - g — 

(1  —  e  cos  u) 


-  e  cos  u  +  v  1  —  e  dL 
e  cos  u)  (1  -  e  cos  u)  dt 


cos  u  + 


3  n(t  -  t  )  sin  u 
2  (1  —  e  cos  u)2 


a  sinu  dL 


(1  —  e  cos  u)' 


-u-  I  cos  u  + 


a  dL 

( 1  —  e  cos  u)  dt 


e  cos  u 


(cos  u  —  e) 
(1  —  e  cos  u) 


(1  -  e  cos  u) 


cos  u  - 


•  2 

sin  u 

( 1  —  e  cos  u) 


d 

(Sal) 

a(cos  u  —  e) 

dL 

d 

/ 9v2  \  a  J 

dt 

\  */ 

j 

(1  —  e  cos  u) 

dt  ’ 

dt 

Ui ;  -  (1 

/ 9v^\  _  1  d2L 

\  9a  )  (1  -  e  cos  u)  ,  2 


(1  —  e  cos  u)' 


2  n(t  -  t  )  (cos  u  -  e) 

-sinu+  j  - ~ — 

( 1  —  e  cos  u) 


,  3n(cosu-e)  dL  1  ,  dL\2 

n  ~  ,3  dt  ,3  '  dt  ' 

(l-ecosu)  (l-ecosu) 

7  n(t  —  t  )  sin  u  ,  ,  , 

x  (cosu-e)  +  - 2 - r  [1+ 

2  (1  —  e  cos  u)  1  (l-ecosu) 

(dA  _  d2L  f  .  3  sinu| 

\9a/  "  (l-ecosu)  ,,2  I cos  u  2  2 

dt  1  ( 1  -  e  cos  u)  - 

,  3n  Ji  -  e2  sinu  dL  \l i  —  e2  ,dL,2 


n  v  1  -  e  sin  u  dL  _  v  1  —  e  ^  dL  ^2 
(l-ecosu)3  ^  (l-ecosu)3  ^ 


..  I  .  ,3 

X  {  sin  u  +  t  Ta _ ~~ _ 


2  (1  —  e  cos  u)  1(1  -  e  cos  u) 


,  .  2 
3e  sin  u 
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t,  (0  - - a^inu  £L  (cqsu+  (cosu^e) 

it"  '  de  /  ( 1  —  e  cos  u)  dt 


£ 

dt 


e  cos 


3 


dt 


ftf) 


(1  -  e  cos  u) 


,2, 

_ a  d  L 

(1  —  e  cos  u)  ^2 


,d L.2 

4  1  dt  ' 


,  i2  ,,  2,  .  3sin2u(l-e2) 

—(cos  u  -  e)  -d-e)+  (T^ecosuT 


e  cos  u 


Jl  - 


\/l  -  eZ 
(i  —  e  cos u) 


cos  u  - 


.  2 

sin  u 


(1  -  e  cos  u) 


( 1  -  e  cos  u) 


4  sin  u  cos  u 


,  dL.2 
3  '  dt 


e  sin  u 
*Jl  -  e2 


-  e2 

(1  -  e  cos  u) 


,  •  3 

3e  sin  u 


(1  -  e  cos  u) 


a(cos  u  —  e)  d^L 
(1  —  e  cos  u)3  dt2 


a  sin  u 


(1  —  e  cos  u) 


d2  /9vi\ 
dt2  V  dI>  ' 

d2  / 3vZ\  _  a  Ji  -  e2  sinu  d2L  a  v/l  -  e2  ,dLx2 

,t2  \~9T /  ,,  ,3  ,  2  .4  dt 

v  '  (1  —  e  cos  u)  dt  (i-ecosu) 


,d_L,2  3e(cos  u  —  e)  , 
4  dt  ‘  (1  -  e  cos  u)  ‘ 


-cos  u  + 


.  .  2 
3e  sin  u 

(1  —  e  cos  u) 


(44) 


(45) 


Let  a  denote  one  of  the  parameters  a,  e,  t.  By  (22)  we  then  have  the  partial  derivatives  of  the 
position,  velocity  and  acceleration  in  Brown's  lunar  orbit  with  respect  to  a,  e,  1  given  by 


d_ 

dt 


SRL 

ii 

2 

cj 

k  95 

k=l 

ii 

"Vl'a 

2 

2 

Cj  -  ( 
c k  dt  \ 

'9vk\ 
^95  ) 

2 

*  2 

dCk  9vk 
dt  95 

k=l 

k=l 

II 

•"Vlia 

2 

z 

cj  d2 

(  9vk> 

)+2 

2  dC^ 
y  Ck  d 

L>  dt  dt 

k=l 

k=l 

+ 

1  d2cJ 

dt2 

„  h 
9v 

95 

j  =  1,2,  3 


(46) 


Next,  we  define 

i  =  inclination  of  Brown's  mean  lunar  orbital  plane  at  time  t  on  the 
mean  ecliptic  of  time  t 

=  ascending  node  of  Brown's  mean  lunar  orbital  plane  on  the  mean 
ecliptic  of  time  tQ  measured  from  the  mean  equinox  of  time  t 
along  the  ecliptic 

=  argument  of  perigee  of  Brown's  mean  lunar  orbit  at  time  tQ 
measured  along  the  orbital  plane  from  the  ascending  node  on 
the  mean  ecliptic  of  time  t  . 
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’**s**r 


Further,  let  i,  £2,  w  be  the  similar  reference  angles  at  the  initial  time  t  of  Brown's  mean  lunar 

12  3  ° 

orbit  relative  to  the  coordinate  system  (y  ,y  ,y  )  referred  to  the  mean  equinox  and  equator  of 
1950.0  as  defined  at  the  beginning  of  this  section.  Preliminary  to  finding  the  partial  derivatives 
of  the  position,  velocity  and  acceleration  in  Brown's  mean  lunar  orbit  with  respect  to  I,  £2,  Z>, 
we  shall  derive  the  expressions  for  the  partial  derivatives  of  iQ,  £2q,  with  respect  to  I,  £2,  w. 
By  (20)  and  the  fact  that  the  matrix  A  is  orthogonal,  we  can  write 

3 

Bjj  =  £  A*  C*  ,  j,  k  =  1,  2, 3  (47) 

1= 1 


where  the  matrix  B  is  given  by  (11)  with  the  angles  i,  £2,  w  replaced  by  iQ,  £2q,  w q,  the  matrix 
C  is  given  by  (11)  with  B^  replaced  by  and  with  the  angles  i,  £2,  w  replaced  by  I,  £2,  Z\  and 
the  matrix  A  is  given  in  Appendix  B  evaluated  at  time  t  .  From  (47)  it  follows  that 


cos  1 


i  =  Yj 


l  l 
C 

o  3^3 

1=1 


0°^  iQ^  180' 


(48) 


sini  sin 
o 


in<.2o=  l  A‘C 


C  »  l„i 
3 

1=1 


—  sin  iQ  cos  £2q  £  A*C* 


1=1 

3 


sin  i  sin  w  =  T  A,C; 
o  o  3  ,1 

1=1 


sini  cosw 


■  l  A3C2 


1=1 


0° <  n  <  360' 
o 


0°^w  <  360‘ 

o 


(49) 


(50) 


Let  a  denote  one  of  the  parameters  i,  £2,  w.  Differentiating  (48)  we  obtain 


9i 

o 

9a 


ac. 


_i —  V  a^  — r- 
sin  i  ^  d  da 


1=1 


Differentiating  the  equations  in  (49),  we  obtain 


(51) 


9£2q  9i 

sini  cos  £2  -r=-  +  cos  i  sin  £2  -r= 

o  c  da  o  o  9a 


-  e  a; 

i=i 


f  9_S 

da 


™  o  9io  v  1  9C3 

sinioSinJ2o  W  -cosi0cos£2o  ^  =  >,  A2  W 

1  =  1 


Multiplying  the  first  of  these  equations  by  cos  £2q  and  the  second  by  sin£2Q  and  adding,  we  see  that 
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(52) 


3  o 
■35 


(cosfi  l  a{  S  4  w)  • 

'a  sin  1  \  o  u  1  on  o  LJ  L  da  ] 

0  \  |=i  1=1  ' 


Differentiating  the  equations  in  (50),  we  obtain 

8"o  t  .  .  8Io  y  Ai  8C( 

smi0  cosWo  g=r  4-  cosio  smw0  =  l  A3  -5=- 

1=1 
3 


3co 


ai. 


<  SC2 


o  o  V  «  1  — 2 

— sini  sin co  -7=-  +  cos  1  coso;  rr=r  =  7,  A,  -r=- 

n  O  uCt  n  n  r>rv  3  rJA' 


o  3a  ^  “3  Da 
1=1 


Multiplying  the  first  of  these  equations  by  coscoQ  and  the  second  by  sinwQ  and  subtracting,  we 
see  that 


9wo  1  / 

-r=-  =  — : — —  COS  CO 

3cv  sini  \  o 


y  A'  !£i_  •  y  A» 

^  A3  3a  Sin'"o  ^  A3  3 a  ) 
1=1  1=1  ' 


(53) 


Finally,  using  expression  (11)  for  the  matrix  C  with  replaced  by  and  the  angles  i,  $1,  co 
replaced  by  i,  SI,  co,  we  see  that 


!Sl _ .2 

30  " 


3ik  _C1 

30 


S 

30 


=  0 


k  =  1,2,3 


(54) 


92i 

Oco- 


Ci 


3C  J 

±1  -_Cj 

3w  1 


3C,J 

=° 


j  =  1,2.3 


(55) 


9C1  3  - 

IT  =  ci  sinJ2 


ocr  3  - 

~W  =~ci  cosa 


3C 


1  .  - 
-TP-  =  sin  co  cos  l 
3i 


9^"2  3  — 

"tt"  =  C,  sin  0 
3i  2 


9C2  3  - 

~W  =-C23cosO 


9C2 

- r=-  =  cos  co  cos  i 
3i 


3C3 

IT 


sin  0  cos i 


3C 


3i 


3C3 

IT 


.3-  =  _ 


cos  0  cos i 


=  —  sm  i 


(56) 


In  (56)  the  angles  I,  0,  co  are  determined  by 

,3 


cos  i  =  Ci 


0”  i  <  180' 


(57) 


sin  i  sin  0  =  C 


-  —  2 
-sin  i  cos  0  =  C3 


0 0  O  <  3  60  1 


(58) 


sin  i  sin  co  =  CJ 


sin  1  cos  co 


- 


0“<:  w  <  360' 


(59) 
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12  3  12  3 

Let  (a^,  cXq,  aj)  denote  the  quantities  (iQ,  ftQ,  wQ)  and  let  (5,3,  a  )  denote  the  quantities 

(T,  ft,  uJ).  Then  by  (22)  the  partial  derivatives  of  the  position,  velocity  and  acceleration  in  Brown's 

mean  lunar  orbit  with  respect  to  the  initial  angles  I,  ft,  w  are 


j  =  1,2,3 
£  =  1,2,3 


(60) 


where  the  3X3  matrix  (da™/dal)  is  determined  from  formulas  (51),  (52)  and  (53)  evaluated  at 
the  initial  time  tQ.  By  (20)  and  (21)  we  have 


12  3 

Let  a  ,  a  ,  a  denote  i,  ft,  w.  Then  by  (6)  we  can  write 


so  that 


+  cv^(t  -  tQ)  +  a£( t  -  tQ)2  +  aj(t  -  tQ)3 


da  ^  da 
o 


£  =  1,  2,  3  ,  k  =  1,  2 


(62) 


2  3 

In  (62)  the  partial  derivatives  with  respect  to  a  =  ft  and  a  =  w  are  given  by  (13)  and  (14),  and 

3 

the  partial  derivatives  with  respect  to  a  =  i  are  given  by 
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L 


9B; 

sin  12 

DB? 

3 

=  -B^  cos  12 

3i 

3i 

=  sin  a;  cos  i 

Ol 

3B2 

T,3 

sin  12 

3B2 

_  3 

!lB2 

9i 

=  B2 

9i 

=  -B2  cos  12 

=  COS  W  COS  1 

01 

By  (12)  and  (62),  we  have 


(63) 


(64) 


£  .  i  c 

where  the  partial  derivatives  of  3B.  / daJ  are  given  by  (13)  through  ( 1 7 )•  with  B.  replaced  by 

f  £  K  R 

0B^/  9a  . 

V.  DETERMINATION  OF  MOTION  AND  PARTIAL  DERIVATIVES 
WITH  RESPECT  TO  INITIAL  CONDITIONS 

1  6 

Let  (x  ,  ....  x  )  denote  the  position  and  velocity  of  the  Moon  relative  to  the  Earth,  and 

I  iiiL.  ^  me 

let  (y  ,  . .  . ,  y^g)  denote  the  position,  velocity  and  acceleration  in  Brown's  mean  lunar  orbit. 
Let 


,k  _  k  _  k 
^me  ~  xme  ^me 


k  =  1, ....  6 


Then  by  (4)  the  (|^  ,  .  .  . ,  |  ^  )  satisfy  the  system  of  equations 

Ail  v*  ill  C 


where 


me  _  k+3 

dt  *  me 


k+3 

a^me 

dt 


=  -yMs 

/M  \ 

Uf 

xme 

3 

k+6  , 
-yme  + 

rme 

+  Hk  + 

{d- 

'  m 

Fk- 

m 

Me  e/ 

some 

t 

.  <+3 
’me 

=  ,  k+3 
’ome 

when  t  =  t. 


k  =  1,2,3 


t  ^  =  (x  k  _  y  k  )i  k  =  l  6 

5ome  'me  yme''t=t  ’  '•••'  • 

o 


(65) 


(66) 


(67) 
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PEP  determines  the  position  and  velocity  of  the  Moon  as  functions  of  time  by  numerically 


integrating  the  differential  equation  system  (66)  for  the  (| 


■ »  Cj  using  relation  (65)  and 


the  fact  that  the  (yme,  ■  -  ,  yme)  are  known  as  functions  of  time.  During  the  integration  the  posi¬ 
tions  of  perturbing  planets  1,  2,  4,  . . . ,  9  in  the  4-k  term  of  (3)  are  determined  from  an  input  mag- 

netic  tape.  The  position  of  the  Earth-Moon  barycenter  relative  to  the  Sun  (x  ,  x  x  )  which 

U  jj  cs  cs’  cs" 

is  needed  in  evaluating  the  B  term  of  (3)  because  of  the  relations 
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k  =  1,  2,  3 


is  determined  either  from  the  perturbing  planet  input  tape  or  from  integrating  the  equations  of 

motion  of  the  Earth-Moon  barycenter  given  in  Ref.  12  along  with  equations  (66). 

Let  9x^e/9^  (j.  k  =  1.  . ,  6)  denote  the  partial  derivatives  of  the  position  and  velocity 

of  the  Moon  with  respect  to  init)al  osculating  elliptic  orbital  elements  (/? 1 , .  .  . ,  /36  ),  and  let 
k  .  —  ■»  ni  m 

3yme/ 9^ni  (j  =  •  ■  • »  k  =  1, ....  9)  denote  the  partial  derivatives  of  the  position,  velocity 

and  acceleration  in  Brown' s  mean  lunar  orbit  with  respect  to  initial  mean  orbital  elements 

<4 . <&>• 


j,  k  =  1,  . . . ,  6 


Then  by  (5)  the  jj  .  (j,  k  =  1,  .  .  . ,  6)  satisfy  the  system  of  equations 
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k  =  1,  2,3 
j  =  1 , .  .  . ,  6  (70) 


±  Fk) 


where 


when  t  =  t 


j.  k  =  1 . 6 


The  9xk  /<>pL  at  the  initial  time  t_  are  determined  from  the  elliptic  orbit  formulas  of  Ref  13 

nit;  III  k  " 

PEP  determines  the  <j.  k  =  1,  .  .  . ,  6)  as  functions  of  time  either  by  assuming  that 

they  are  equal  to  the  3y  k  / 90  ^  (j,  k  =  1,  .  .  . ,  6)  or  by  numerically  integrating  the  differential 
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k  k 

equation  system  (70)  for  the  n  (j,  k  =  1, . . . ,  6)  using  relation  (69)  and  the  fact  that  the  9y  / 

.  me 

dp J  (j  =  1, ....  6,  k  =  i,  . . . ,  9)  are  known  as  functions  of  time.  PEP  has  these  two  options, 

m  k  /  i  k  /  -  i 

because  it  might  be  sufficiently  accurate  to  assume  that  3x  /3 p*  =  3y  /3 pJ  in  the  least- 

°  me  m  me  m 

squares  process  of  fitting  to  observations.  As  mentioned  earlier,  not  having  to  follow  the  exact 
procedure  of  numerically  integrating  equations  (70)  would  save  a  great  deal  of  computer  time. 

We  intend  to  adjust  the  initial  osculating  elliptic  orbital  elements  rather  than  the  initial 
position  and  velocity  in  the  least-squares  process  of  fitting  ephemerides  to  radar  and  optical 
observations,  because  the  method  of  integration  we  use  makes  it  easy  to  calculate  the  partial 
derivatives  of  position  and  velocity  with  respect  to  the  former  parameters  and  because  the  former 
parameters  should  be  less  correlated  than  the  latter.  The  elliptic  orbital  elements  (/?*,  . .  . ,  p^) 
we  have  used  in  this  report  and  in  Ref.  i  are: 

a  =  semimajor  axis 
e  -  eccentricity  (0  ^  e  <  1) 
i  =  inclination  (0“<:  i<:  180°) 

&  =  longitude  of  ascending  node  (0°^  £1  <  360°) 

w  =  argument  of  perigee,  measured  along  the  orbital  plane  from  the 
ascending  node  (0  4  w  <  360  °) 

i  =  initial  mean  anomaly 


Some  other  choice  of  elliptic  orbital  elements  could  be  less  correlated  l  elative  to  the  radar  and 
optical  data  than  this  choice.  For  instance,  this  would  very  likely  be  the  case  if  we  replaced 
the  angles  (£2,  w,fQ)  by  the  angles 


£2  =  £2 
5  =  £2  +  w 


M  =  £2  +  w  +  £ 

o 


(72) 


since  the  angles  in  (72)  are  those  generally  used  in  celestial  mechanics  and  were  probably  arrived 
at  from  the  results  of  experience  with  optical  data.  The  formulas  in  this -report  and  in  Ref.  1  are 
stated  in  terms  of  the  angles  (£2,  w,  t  ).  The  partial  derivatives  with  respect  to  (?2,  Z,  M)  can  be 
determined  from  these  results  because  of  the  relations 
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(73) 
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APPENDIX  A 

EXPANSION  FOR  PRECESSION  MATRTX 
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Let  iy*.y2.y3)  be  a  coordinate  system  referred  to  the  mean  equinox  and  equator  of  1950.0 

12  3 

(J.  E.  D.  2433282.423),  and  let  !x  ,  x  .  x  )  be  a  coordinate  system  referred  to  the  mean  equinox 
and  equator  of  date.  Then  the  relation  between  them  is  given  by 


yJ  =  V  Pjxk 
J  '  k 


J  =  1.2.3 


J  V  ,,  k  k 
xJ  =  P  y 

k  .=  1  J 

where  the  orthogonal  matrix  (P^)  is  the  precession  matrix. 

k  14 

To  give  the  established  expression  for  the  precession  matrix,  we  first  define  the  angles 

r  -  2304 V948T  +  0V302T2  +  0VO179T3 
o 

z  =  2304'.’948T  +  1V093T2  +  0V0192T3  '  <A_2 

G  =  2004 '.'255T  -  0V426T2  -  070416T3  . 

where  T  is  measured  in  tropical  centuries  of  36524.21988  ephemeris  days  from  the  epoch 

1950.0  (J.  E.  D.  2433282.423)  to  the  instant  of  interest.  Then  the  precession  matrix  at  this  in- 

.  •  15 
stam  is 

P  }  =  cos  f  cos  6  cos  z  —  sin  K  sin  z 
1  *  o  o 

P?=  —  sin?  cos  G  cos  z  —  cos  t  sinz 

1  •  o  '  o 

P  3  =  —sin  Q  cos  z 

P  3  =  cos  f  cos  G  sin  z  +  sin  £  cos  z 

2  ‘  o  o 

p.?  =  -sint  cos  G  sinz  +  cosf„  cos  z  }  .  (A-: 


P^  =  —sin  G  sin  z 


P,  -  cos  f  sin  G 
3  o 


P  =  —  sin  £  sin  0 
3  o 


P3  =  cos  G 


Let  r  denote  the  time  from  the  epoch  1950.0  (J.  E.  D.  2433282.423)  in  units  of  10,000 
ephemeris  days.  Then  by  Taylor's  theorem  we  have 


PRECEDING  PAGE  BLANK 
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Treating  the  coefficients  in  (A-2)  as  exact,  some  simple  calculations  show  that  the  terms  up  to 
the  fifth  power  in  the  Taylor  expansions  (A-4)  are: 

P,1  =  1.0  -  2.22603398052517  X  10_5t2  -  2.6903385325366  X  10_9r3 
1 

+  8.191221606878  X  10_11t4  +  1.79948222850  X  10-14t5 
P  2  =  -6.119064710033514  X  10_3t  -  5.06975739290688  X  10_7r2 

+  4.5321716219079  X  10_8r3  +  8.619581795926  X  10_12r4 

-  1.02943658327  X  10-13r5 

P3  =  -2.660399722772102  X  10-3t  +  1.54818397804898  X  10"7r2 

+  1.9729201591810  X  10-8r3  +  1.960730253191  X  10-12t4 
-4.39298354075  X  10~14r5 

P^  =  6.119064710033514  X  10“3t  +  5.06975739290688  X  10-7-2 

-  4.5321716219079  X  10-8r3  -  9.636891635856  X  10-12/1 
+  1.02604298897  X  10~13r5 

.  (A-5) 

P 2  =  1.0  -  1.87214764627888  X  10~5r2  -  3.1022173551368  X  10_9t3 

+  6.882478825535  X  lO"1^4  +  1.91215207447  X  10_14t5 
P 3  =  -  8.13957902909886  X  10_6r2  -  5.8309700675934  X  10"10r3 

+  2.994360606802  X  10-11-4  +  5.71739459043  X  10_15r5 
P^  =  2.660399722772102  X  10-3r  -  1.54818397804898  X  10"?t2 

-  1.9729201591810  X  10"8r3  -f  3.791379581151  X  10"13t4 
+  4. 50404085077  X  10-14t5 

P32  =  -8.13957902909886  X  10-6t2  4  1.8168268497009  X  10-10t3 

+  3.024323052660  X  10-11r4  +  2.58550054981  X  10-17t5 
•P3  =  1.0  -  3.53886334246294  X  10_6r2  +  4.1187882260017  X  10_1°t3 

+  1.308742781343  X  10-11t4  -  1.12669845971  X  10-15t5 
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APPENDIX  B 

EXPANSION  FOR  MATRIX  A  AND  ITS  DERIVATIVES 


The  mean  obliquity  of  the  ecliptic  is 


16 


e  =  ]3°27l08l.'26  — 46'.,845T-0'.'0059T2  +  0V00181T3 
o 


(B-l) 


where  T  is  measured  in  Julian  centuries  of  36525  ephemeris  days  from  the  epoch  1900  January 

0.5  E.  T.  -  J.  E.  D.  2415020.0  to  the  instant  of  interest.  The  relation  between  a  coordinate  sys- 
12  3 

tern  (x  ,  x  ,  x  )  referred  to  the  mean  equinox  and  equator  of  date  and  a  coordinate  system 

123  *17 

(w  ,  w  ,  w  )  referred  to  the  mean  equinox  and  ecliptic  of  date  with  the  same  origin  is  ‘ 


1  1 

x  =  w 

2  2  3  . 

x  =  w  cos  £  —  w  sin  c 
o  o 

3  2  .  .3 

X  =W  SI. If  +W  COS  £ 

o  o 


(B-2) 


Let  t  denote  the  time  from  ‘he  epoch  1950.0  (J.  E.  D.  2433282.423)  in  units  of  10,000  ephemeris 
days.  By  treating  the  coefficients  in  (B-l)  ar.  exact,  Taylor  expansions  similar  to  (A-4)  gne 


sin  c  =  0.3978811865927521  -  5.70513893192403  X  10_4t 
o 

-  1.8312087506169  X  10-9t2  t  1.652267540061  X  10_1°73 

•r  4.45783951328  X  10“15r4  -  2.36469209  X  10_19r5 

cos  £  =  0.9174369522509674  +  2.47424898500217  X  10"5r 
o 

-  1.31335717*0992  X  10-9r2  -  7.173527734648  X  10_1173 
+  1.02732897621  X  10-14t4  i  3.29806267  X  10-1975 


(B-3) 


Combining  transformations  (A-l)  and  (B-2)  gives  transformation  (18),  where 


Aj  =  Pi  cos  e  +  Pi  sin  e  ■ 

2  2  o  3  o 

Ai  =  Pi  cos  £  -  Pi  sin  e 

3  3  o  2  o 


j  =  1,2,3 


(B-4) 


The  Taylor  expansions  for  the  A^  (j  =  1,  2,  3)  are  given  by  the  first  three  equations  in  (A-5).  To 
determine  the  Taylor  expansions  for  the  A^,  (j  -  1,  2  3),  we  multiply  and  add  the  expansions 
given  in  (A-5)  and  (B-3)  as  indicated  in  (B-4): 
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A^  =  6.672379076707188  X  10-3r  +  4.03140345445988  X  10_7r2 

-  4.9421227156709  X  10“8r3  -  8.685948302421  X  10_12r4 
+  1.11902061670  X  10"13t5 

A2  =  0.9174369522509674  +  2.47424898500217  X  10"5t 

-  2.04156730272966  X  10_5r2  -  2.8443776398581  X  10-9r3 
+  7.513826113715  X  10-11t4  +  1.75327407517  X  10_1475 

A3  =  0.3978811865927521  -  5.70513893192403  X  10"57 

-  8.87742893170170  X  10_672  -  2.0534553328574  X  10-1°r3 
+  3.266231487901  X  10_11t4  +  4.79023554313  X  10-15t5 

A^  =  6.08828576338671 X  10_67  4  7.11738284222153  X  10-8r2 

-  3.4836043645777  X  10_11t3  -  9.238939614350  X  10_1474 

-  1.72691125954  X  10~1675 

A2  =  -0.3978811865927521  +  5.70513893192403  X  10_5t 

-  1.67960985290526  X  10-8r2  -  3.37101 16720406  X  10-11t3 
+  1.616276978960  X  10'l3r4  +  7.61970380671  X  10_16r5 

A3  =  0.9174369522509674  +  2.47424898500217  X  I0'5T 

-  9.41199405263501  X  10_9r2  -  1.3793679110716  X  10_11t3 
+  6.983259897028  X  10"14r4  +  3.21079987643  X  10"16j5 


(B-5) 


The  expansion  for  the  matrix  A  given  in  (B-5)  and  the  first  three  equations  of  (A-5)  has  13 
decimal  place  accuracy  30  years  away  from  the  epoch  1950.0  and  9  decimal  place  accuracy  300 
years  away  from  the  epoch  1950.0.  Actually,  we  treat  expansions  (A-5)  and  (B-5)  as  defining 
the  matrix  A  so  that  the  formulas  for  the  mean  lunar  orbit  involving  the  matrix  A  can  be  r  egarded 
as  exact,  except  foi  those  which  assume  that  the  inverse  of  A  is  equal  to  its  transpose.  But 
this  is  only  done  in  computing  the  matrix  [da  da^)  in  (60),  so  that  only  those  forrm  1  s  involving 
partial  derivatives  with  respect  to  initial  mean  orbital  elements  (i,  S2 ,  u>)  are  affected.  Because 
of  the  use  which  is  to  be  made  of  these  quantities,  any  possible  loss  in  accuracy  due  to  the 
assumption  that  the  inverse  of  A  is  equal  to  its  transpose  is  unimportant. 

Let  t  denote  time  measured  in  days.  For  a  function  f  defined  by 

5 

f  =  y.  f  rn 

u  n 
n=0 
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u 


we  have 


2  2 

Thus,  the  coefficients  in  the  expansions  for  dA/dt  and  d  A/dt  are  easily  derived  from  the 
expansion  for  the  matrix  A  given  in  (B-5)  and  the  first  three  equations  of  (A-5). 


23 


PRECEDING  PAGE  BLANK 


APPENDIX  C 

DIFFERENCE  BETWEEN  TRUE  LUNAR  ORBIT  AND  MEAN  LUNAR  ORBIT 

The  graphs  of  Fig.  C-l(a-c)  represent  the  differences  between  the  x,  y  and  z  components 
and  the  radius  vectors  in  the  true  lunar  orbit  and  in  Brown's  mean  lunar  orbit  during  the  •'50- 
dav  period  from  9  January  1967  to  3  April  1968.  The  distance  unit  is  earth  radii  and  the  coor¬ 
dinate  system  is  referred  to  the  mean  equinox  and  equator  of  1950.0.  The  true  lunar  orbit  coor- 

1 8 

dinates  were  taken  from  the  Jet  Propulsion  Laboratory  Ephemeris  Tapes;  the  mean  lunar  orbit 
coordinates  were  evaluated  by  using  the  formulas  in  tnis  report. 

From  these  graphs  it  is  clear  that  the  mean  lunar  orbit  does  indeed  follow  the  true  orbit  on 
the  average  as  the  Moon  moves  at  a  distance  of  60  earth  radii  from  the  Earth.  However,  there 
are  considerable  oscillations  about  the  mean  orbit  (mostly  due  to  the  Sun),  so  that  numerically 
integrating  the  equations  for  the  difference  between  the  true  and  mean  lunar  orbits  represents 
only  a  1  ]-  order-of-magnitude  improvement  over  numerically  integrating  the  original  equations 
of  motion.  Althougn  tnis  is  a  considerable  saving,  it  is  not  as  large  as  one  might  hope. 

Brown's  lunar  theory  represents  the  motion  of  the  Moon  by  the  mean  lunar  orbit  plus  over 
1650  trigonometric  terms.  Adding  a  few  of  the  larger  trigonometric  terms  to  the  mean  lunar 
orbit  would  yield  an  orbit  which  gives  a  further  saving  in  representing  the  true  orbit  of  the  Moon. 
(From  the  graphs  it  would  appear  that  the  main  terms  of  this  addition  would  have  periods  of  about 
14  days  and  somewhat  less  than  a  year.)  At  some  future  time  we  shall  alter  the  mean  lunar  orbit 
subroutines  in  PEP  in  this  manner;  for  the  present  they  utilize  the  formulas  in  this  report. 
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